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We develop a flux-conservative formalism for a Newtonian multi-fluid system, including dissipation 
and entrainment (i.e. allowing the momentum of one fluid to be a linear combination of the velocities 
of all fluids). Maximum use is made of mass, energy, and linear and angular momentum conservation 
to specify the equations of motion. Also used extensively are insights gleaned from a convective 
variational action principle, key being the distinction between each velocity and its canonically 
conjugate momentum (which is modified because of entrainment). Dissipation is incorporated 
to second order in the "thermodynamic forces" via the approach pioneered by Onsager, which 
makes it transparent how to guarantee the law of increase of entropy. An immediate goal of 
the investigation is to understand better the number, and form, of independent dissipation terms 
required for a consistent set of equations of motion in the multi-fluid context. A signiflcant, but 
seemingly innocuous detail, is that one must be careful to isolate "forces" that can be written as total 
gradients, otherwise errors can be made in relating the net internal force to the net externally applied 
force. Our long-range aim is to provide a formalism that can be used to model dynamical multi- 
fluid systems both perturbatively and via fully nonlinear 3D numerical evolutions. To elucidate the 
formalism we consider the standard model for a heat-conducting, superfluid neutron star, which is 
believed to be dominated by superfluid neutrons, superconducting protons, and a highly degenerate, 
ultra-relativistic gas of normal fluid electrons. We determine that in this case there are, in principle, 
19 dissipation coefficients in the final set of equations. A final reduction of the system is made by 
neglecting heat conduction. This leads to an extension of the standard two-fluid model for neutron 
star cores, which has been used in a number of previous applications, and illustrates how mutual 
friction is represented in our formalism. 

PACS numbers: 



I. INTRODUCTION 

The chemistry community has great experience in dealing with multi-fluid systems. Perhaps indicative of the 
cultural divide between research areas is the development of a multi-fluid literature for superfluid neutron stars which 
largely ignores the chemistry successes. Of particular note is the Onsager formulation of dissipative, multi-fluid 
systems [1]. On the other hand, action-based derivations of the equations of motion for multi-fluid systems [2-5] 
have, as far as we can tell, not crossed over in the other direction. In this work, we attempt to span the divide by 
utilizing fully mass, energy, and linear and angular momentum conservation, the Onsager formulation, and action- 
based constructs to build a general dissipative model for multi-fluid systems. We facilitate our use of the conservation 
laws by borrowing from the engineering community the idea of "control" volumes. Even though we restrict the 
discussion to the Newtonian regime, we do not hesitate to use ideas that have been developed for general relativity. 
In fact, the analysis is often much simplifled if the true covariant nature of the physics is retained. 

Although much of our analysis is general, and can be used in any multi-constituent fluid context — for example, 
superfluid Helium and bubbly liquid/gas mixtures (see [6] for a general review of two-phase flow, and a series of 
papers by Geurst [7-9] for interesting discussions), our immediate goal is to produce a formalism that can be applied 
to Newtonian superfluid neutron stars, in order to model their detailed dynamics. The importance of this target 
is enhanced by the fact that the first generation of large-scale, ground-based, gravitational-wave detectors are now 
operating near their design sensitivity [10]. Mature neutron stars can radiate gravitationally in a number of ways, 
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many of which depend crucially on dissipation in the superfluid interior. An obvious example is the potentially strong 
damping of gravitational-wave driven mode-instabilities (eg. of the so-called r-modes) due to either superfluid mutual 
friction [11, 12] or a viscous boundary layer at the neutron star core-crust interface [13, 14]. Other astrophysical 
contexts requiring a detailed model for superfluid dissipation include the pulsar glitches [15], the standard model for 
which is based on the transfer of angular momentum between components rotating relative to each other, and possible 
neutron star free precession [16], the damping of which depends crucially on the interior viscosity. 

In the standard superfluid neutron star model, the neutrons in the inner crust and the outer core are believed to be 
superfluid, and the protons in the outer core are thought to be superconducting. Because of overall charge neutrality, 
the neutrons and protons will coexist with a normal, but highly degenerate fluid of electrons (and possibly muons in 
the more dense regions) . At the core temperature expected in a mature neutron star (likely four orders of magnitude 
below the relevant Fermi temperatures), the superfluid and superconducting nature, respectively, of the neutrons and 
protons implies that they function as dynamically independent, interpenetrating fluids. The conventional wisdom is 
that viscous damping of the oscillations of these objects is a result of electron-electron scattering. However, this begs 
the following question: how is it that damping of the electron motion on the one hand leads to damping, on the other 
hand, of the superfluid neutrons, which make up the bulk of the matter in the star? Despite these kinds of problems 
having been studied for decades [17, 18], the answer is not at all obvious from the extant literature (although we 
have recently clarifled some aspects in [19, 20]). Indeed, the conventional wisdom seems to have been guided by the 
classic discussion of Cutler and Lindblom [21], who noted that "neutron star matter becomes more viscous in the 
superfluid state than it was in the normal state." However, while they included the effect of superfluidity in their 
determination of the viscosity coefficients, they did not consider the effect on the actual fluid dynamics. Obviously 
aware of this fact, they remark that "The superfluidity of neutron star matter, will, of course, have drastic effects both 
on the dissipation coefficients and on the dynamics of this material." The implications of the first statement on the 
possibility of detectable gravitational radiation from neutron stars are potentially far-ranging, leading most directly 
to the conclusion that gravitational- wave driven instabilities are ineffective in mature neutron stars [11], provided, of 
course, that superfluid dynamics can be ignored. However, subsequent work on the role of superfluidity in neutron 
stars reveals that the effects on the dynamics can be significant, and sometimes unanticipated (see [22-25] for recent 
discussions). 

To date, most studies of superfluid neutron stars have drawn upon Landau's original two-fluid model for He^ 
[26]. In the case of neutron stars, both neutrons and protons are expected to form Cooper pairs, which then form 
condensates which at the fluid level can be described by the standard equations. The main distinction in the neutron 
star case is that the two fluids no longer represent the superfluid and normal parts of a single particle species. Instead, 
the two degrees of freedom describe the neutrons and a conglomerate of all charged components which are expected 
to be electromagnetically coupled on a short timescale. Following the traditional route, Mendell [18, 27] extended 
the non-dissipative zero-temperature equations to account for the main dissipation mechanisms. The result is a set 
of equations which include the mutual friction coupling due to electrons scattering off of rotational vortices in the 
condensates (see also [20]). These equations were later used to study the damping of unstable modes of oscillation in 
a spinning neutron star. The conclusion drawn from this work is that mutual friction is a key dissipation mechanism 
which may, in fact, suppress rotational instabilities entirely [12, 28]. 

As the work of Mendell shows, the mathematical description of the problem becomes considerably more complicated 
once the so-called entrainment effect [29] , which accounts for the fact that the flow of one fluid may impart momentum 
in the other fluid, is included. This complexity is partly due to a common confusion concerning Landau's description 
of superfluid hydrodynamics. At the heart of the problem lies a failure to distinguish clearly between transport 
velocities and momenta. This distinction is made clear in recent variational formulations of the problem [2-5] , work 
which draws heavily on the fully covariant relativistic description of Carter and collaborators [30-36], and which 
incorporates entrainment in a natural way. This leads to a mathematical description which is in many ways less 
complex than the standard one. 

Our present analysis differs somewhat from the recent variational derivations [2-5] by emphasizing the conservations 
laws for mass, linear and angular momentum, and energy. We flnd this approach appealing because it stays close 
to our physical intuition. Furthermore, we believe that our discussion may serve as a useful "introduction" to 
the more technically involved work described in [2-5]. Of course, we rely heavily on the variational methods for 
rigorous mathematical support of our results in the non-dissipative case. Until very recently there had been no 
serious discussion of the dissipativc problem within the variational framework. The discussion by Carter and Chamel 
[3-5] provided a useful first step by translating (and extending somewhat) the covariant relativistic model devised by 
Carter [30] and Carter and Khalatnikov [37], yet the discussion remains somewhat abstract. A key aim of the present 
discussion is to make contact with applications and prepare the ground for building realistic models of astrophysical 
supcrfluids. 

The present work has two main aims. First of all, we want to consider the multi-fluid dissipation problem in 
"complete" generality. In doing this we rely on the key conservation laws and ask what the most general form for the 
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dissipation and interaction forces may be. The dissipativc terms are then constrained by ensuring that the second law 
of thermodynamics is satisfied. To make progress in this direction we make use of Onsager's celebrated symmetry 
principle [1]. Our derivation follows the guidelines set out in work on superfluid He^ (see [38] and [39] for nice 
descriptions) and is basically an extension of the general analysis championed, for instance, by Landau and Lifshitz 
[26] for one-fluid systems. We extend previous work in two important ways: we include the entrainment effect in 
a (hopefully) clear way, and we also allow for the presence of mutual interaction terms which have traditionally not 
been considered (eg. in [3-5]). The latter turns out to be crucial if we want to be able to represent the mutual friction 
forces which arise due to the presence of vortices. Up to this point, i.e. in Sees. II and III, the discussion is general 
and should be relevant for any multi-fluid system. Our second objective is met in Sec. IV, where we specify the 
equations to the conditions that should prevail in the outer core of a mature neutron star. We consider a mixture of 
superfluid neutrons, superconducting protons, normal electrons and entropy. This model is instructive in a number 
of ways. We show how the dynamical degrees of freedom are reduced in a natural way if we assume that the electrons 
and protons are coupled electromagnetically. (It should be noted that we do not account for detailed electromagnetic 
effects in this analysis [40].) We also deduce that this system in general requires the specification of no less than 
19 distinct dissipation coefficients. This discussion is followed, and the paper concluded, by a final reduction to the 
dissipative extension of the standard two-fluid model which has been used in recent work on neutron star dynamics. 
This model, which is obtained by neglecting heat conduction, provides the simplest system that should be considered 
in "realistic" applications. We discuss the interpretation of the various dissipative terms, some of which appear to 
be new. In particular, we show how the standard form for the mutual friction force [20] is represented within our 
formalism. 



II. CONSERVATION LAWS FOR MULTI-FLUID SYSTEMS 

Any general fluid formalism must respect the three laws concerning 

• Conservation of Mass 

• Conservation of Linear and Angular Momentum 

• Conservation of Energy 

Although entropy is not conserved in general, the law of increase of entropy is also a vital property of any fluid. In 
imposing these conservation laws for a multi-fluid system, we consider some "control volume" V with boundary dV, 
letting rji denote the unit normal form to dV. The control volume is to be contrasted with the notion of a fluid 
element, which is a region of the fluid small enough to be considered "infinitesimal" with respect to the whole fluid 
and yet contain enough particles that a (local) thermodynamic treatment is warranted. 

We consider a general multi-fluid system, the dynamical description of which is based on determining the space and 
time dependence of each particle number density and the corresponding particle number density current n^. These 
are natural variables to use in describing the system. It is assumed that the different fluids are interpenetrating, 
i.e. the different constituents do not have to be separated by interfaces. We let x, y, z be particle constituent indices 
that range over the whole set of fluids in the system. Repeated constituent indices will not satisfy a summation 
convention, but vector indices i,j,k do satisfy such a convention (as usual). Ostensibly, all formulas will be written 
in a general coordinate basis, with gij = gji representing the (flat-space) metric, except for certain integrals written 
below, where a Cartesian basis is implicit (since tensor components are being integrated). The metric is used to raise 
and lower indices, e.g. for a vector Vi = gijV^ and = g^-'vj, where the inverse metric g*-' satisfies g^^gui = 5^ i and 5^ i 
is the Kronecker-delta tensor. When integrating over volumes, the measure is given by dV = y/gd^x where g is the 
determinant of the metric. Finally, all spatial derivatives are given by the covariant derivative Vj that is compatible 
with the metric (i.e. Vigju = 0). 

A. Conservation of Mass 

Let denote the "particle" mass of the x*'^-constituent, so that = m^rix is the mass density of that same 
constituent. The total (local) mass density of the system, p, is thus 



X X 



(1) 
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The total mass-density current p* is likewise given by 

p^ = 5^mX- (2) 

X 

When p is integrated over the control volume V, we obviously determine the mass of fluid M contained in V, i.e. 

M = [ pdV. (3) 



When is integrated over dV, we have a representation of the amount of mass that leaves (or enters, depending on 
the situation) V. Overall conservation of mass thus means 

^ [ pdV = - [ p^VidA = - I Vip' dV . (4) 
cir Jv Jav Jv 

This implies the local condition 

dtp + Vip' = . (5) 
If Fx represents the particle number creation rate per unit volume of the x^'^-constituent, then we have 

dtn^ + V,< = Fx . (6) 
Multiplying this by m^, summing over all x, and using the conservation of mass, Eq. (5), one immediately sees that 

J2m^T^ = 0. (7) 

X 

Worth noting at this point is that we will consider the entropy to be a fluid, the main distinction being that its mass 
is zero. Its "number" density is simply the entropy per unit volume, and it is assumed to satisfy a creation rule like 
Eq. (6). 



B. Conservation of Linear Momentum 



Let represent the (local) linear momentum per fluid element (which will later be understood to be canonically 
conjugate to n* ) of the x*^-constituent. The total (local) linear momentum of the fluid Pi will then obviously be 

given by 

Pi = J2Pi- (8) 

X 

Likewise let tt^ denote the total linear momentum density, so that 

7ri = ^nxP^ = ^< . (9) 

X X 

The total (global) linear momentum Pj of the fluid contained in the control volume will be given by 

Pi= f TTidV . (10) 

Jv 

Finally, let T^i represent the t*'^-component of the total linear momentum that flows (per unit time and unit area) 
in a direction perpendicular to the j^^-dircction, and let fi be the net external force density acting on the control 
volume. Overall conservation of linear momentum in V is then given by 

which implies that 



(12) 
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Here we could have added a term that corresponds to the net internal forces. However, the weak law of action and 
reaction implies that this contribution is zero [41]. There are also the well-known difficulties for action and reaction 
when charged particles are present (of. [41]). 

The individual fluid components will not satisfy exactly the same conservation law, since momentum and energy 
can be exchanged between them. However, it is worthwhile to assume that the x*'^-constituent satisfies an equation 
of the form 

dtTTf + VjT-^i = 7f , (13) 

noting that this is still completely general, as nothing has been proposed about either T^^i or ff. That being said, 
one should not be too surprised by such a choice, since it represents the limiting form familiar from single-fluid studies. 
Consistency with overall linear momentum conservation now demands that 

Although it may seem natural to suppose that the left- and right-hand sides should vanish independently, we will 
establish below that this is not correct. 



C. Conservation of Energy 
Let U denote the total energy per unit volume in V, and U be the total energy of the control volume, i.e. 



= / UdV . (15) 
Jv 



U= / UdV 
'v 

Also let denote the flow of energy (per unit time and unit area) perpendicular to the z*''-direction, and e'^''* the 
energy creation rate per unit volume due to external sources. Conservation of energy in V is then naturally expressed 
as 



dU _ d_ 
d^ ~ dt 



f UdV = - [ Q% dA+ [ e'^'^'dV = - f (ViQ' - e<^'^*) dV , (16) 

Jv JdV Jv Jv 



or 



dtU + ViQ' = e'^"* . (17) 

Here it should be noted that we must specify a relation U = U{n^, n^) to close the system of equations. That is, we 
need to provide an equation of state for the matter. 



D. Conservation of Angular Momentum 

In order to get an appropriate definition of the total angular momentum L' in the control volume, let us first 
reconsider the total linear momentum in the context of Newton's Second Law. We need to determine those 
internal and external contributions that act within the fluid to change Pj, so as to properly define the net torque. 

Referring back to Eq. (10), we can use Pj in Newton's Second Law to determine the net "force" Fj acting on the 
control volume: 

= ^ = ^ dtTTi dV = {fi - VjT^i) dV . (18) 

This suggests that the appropriate definition of the total torque T' acting on the control volume (in Cartesian 

coordinates) is 

= ^e'^" {xj [fk - ViT^k] - Xk [fj - ViT^j] ) dV . (19) 
To see that this definition is consistent, we consider the time rate of change of L* and show that it is equal to T\ 
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The total angular momentum of the control volume is defined to be 

V = ^e'^^ (x.-TTfe - Xk-Ki) dV , (20) 

and its time rate of change is equal to 

dL' 1 f 

"dT " 2^''^ j^i'^A'^k-Xkdtnj) AV 

= {xj [h - ViT\] - xk [fj - ViT'j] ) dV 

= T , (21) 

which is the desired result. Of course, Eq. (21) implies that the total angular momentum contained within the 
control volume is conserved when the net torque vanishes. We can use this to our advantage now, by showing that 
conservation results only if the symmetry = Tji holds. 

When a system is completely isolated, and the control volume is such that it completely contains the system and 
has a boundary taken well outside the system, we should find that L' is conserved, i.e. constant in time, if in addition 
the net force fi acting is zero. Eq. (21) can be rewritten as 



^ = l^''' [j^ i^jf" - ""kfj) (Tjk - Tkj) dV - j^^ {x,T\ - x.T'j) m d^l . 



(22) 



The first term on the right vanishes because the net force fi is zero. The last term vanishes because wc have taken 
the boundary of the control volume well outside the system. Hence, if angular momentum is to be conserved, the 
middle term must also vanish, which implies Tij = Tji. 

Because of the inherent linearity in our definitions when summing over the constituent indices, we immediately see 
that an appropriate definition for the control volume force acting on the x*'^-constituent is 

F- = / if- - ^jT-h) dV , (23) 
Jv 

while for the control volume angular momentum we have 

K=le^^' l^{xjnl-Xkn^)dV , (24) 

where the corresponding torque is 

= r''' (^^- Ul - VzT-'fe] - Xk [7] - V,T-',]) dF . 



(25) 



This naturally implies 



di4 

dt 



= l^^'" {x,ft - Xkf^) dV-J^ {T-jk - T'^kj) dV - j^^ {xjT-\ - XkT-'j) m d^j 



(26) 



At this point we make an important observation: In contrast to the comments following Eq. (22), we argue that it is 
a mistake to infer that T^tj = T^n ■ This assumption would be too restrictive. In particular, we will show later that 
it would not allow us to represent the standard form for the vortex-mediated mutual friction in a superfluid system 
[20]. 



E. Lagrange's Generalized Action Principle and Equations of Motion 

Further progress can be made in specifying the unknown quantities U, Q*, /f , T^^ i and i, assuming that fi and 
e'*'^* are prescribed. Consider that there is a Lagrangian density C whose variation leads to the multi-fluid equations 
in the fully conservative, non-dissipative, case. It consists of the total kinetic energy density minus the "potential" 
energy density, which we take to be the total internal energy density £(nx, n]^, 

^ = 11 l^^'aaKK/n^ - S{n., <) . (27) 
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Minimizing the action associated with C in the usual way will necessarily lead to the conservative Euler-Lagrange 
equations. To incorporate dissipation, we use Lagrange's generalized approach (which is championed by Carter, see 
for example [30]), which is to obtain the Euler-Lagrange equations in the usual way, but instead of equating them 
to zero we put them equal to an appropriate sum of generalized forces that may include dissipation. The net effect 
is that the Lagrangian serves to define the momenta, whose identification is perhaps the most crucial element of the 
multi-fluid construction that fully incorporates entrainment. 

The first step towards identifying the unknown quantities is to write the total energy density U (to be distinguished 
from the total internal energy density £) in terms of the fundamental variables of the action. Recall that a Hamiltonian 
density H can be generated as a Legendre transformation on C, i.e. 



(28) 



The momenta are defined via 



dL 
dni 



d£ 
dni 



(29) 



where = n^^/n^ is the transport velocity of the x^'^-constituent. The simplest, and most natural, way to obtain 
the total energy density in terms of the fundamental variables is to setU =H. 

Accepting that, we now take a time-derivative of H, and use Eqs. (6) and (13), to find 



+ 



H +-m v^- v'^Pi 



Fx , 



(30) 



where = gijvl-vl and 



d£ 
Srix 



(31) 



Comparing with Eq. (17) we can identify 



(32) 



1 



+ 



(33) 



It is useful at this point to note that there is a more complete Legendre transformation on £ that can be made, by 
using the functions pQ defined as 



dc 



= - M 



(34) 



This leads to the new function given by 



ii.^-^[n^p-,+nipt]-cj . 



(35) 



We will now establish that enters the total momentum conservation equation as a generalized pressure. To do this 
we need to analyze more closely the stress-tensor T^ j. 
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Widely discussed in the general relativity literature (for instance, Ref. [42]), is the fact that the stress-energy- 
momentum tensor can be obtained from a variation of the Lagrangian with respect to the metric. It is natural 
that metric variations are associated with stresses, since stress leads to deformations and the metric is what gives 
distances and volumes. In Appendix A we provide a non-rigorous account (using energy conservation arguments) of 
how to obtain the conservative piece of T^i, to be denoted C^i, by varying C with respect to Qij (see [2] for a rigorous 
derivation). Our variation is to be understood as a displacc;ment, in a small amount of time, of each fluid element 
within a fixed control volume, keeping the individual particle numbers and their associated velocities fixed. 

The discussion in Appendix A prompts us to consider a less general form of the internal energy, i.e. 

£ = £{n^,wl^) , (36) 

where 

<y =K-Vy , t«xy = gijKyKy ■ (37) 

This is manifestly Galilean invariant and locally isotropic. It is worth pointing out that this is not the most general 
form for the energy that we can allow. For example, it would be straightforward to include an explicit dependence 
on the vorticity, eg. in terms of the "vorticity" vector W^s (defined below by Eq. (60)) we could allow the energy to 
depend on W^g = gijW^^Wl^. This would be analogous to the superfluid Helium model discussed by Bekarcvich and 
Khalatnikov [43] where the circulation around a vortex contributes to the energy budget, and would lead to terms 
that could be interpreted as the "vortex tension" . Although one can think of other such possibilities we prefer to 
work with Eq. (36) here. Once the complete framework has been developed for this case, it should be relatively easy 
to extend it to include more general cases. 

From Appendix A we then find that i is given simply as 

Ch = ^5h + Y.Kpt . (38) 

X 

Thus, ^ behaves as the total pressure. We are naturally led to write the full stress-tensor as a sum of C^i and a 
dissipation piece D^i, i.e. 

= Ch + . (39) 

At this point it is important to note that, since the pressure ^ is not generally separable there is in general no way to 
write i as just a sum of the T^^ i. In other words, there is in general no useful notion of partial pressures, a point 
which is examined in more detail in Appendix B. Basically, it is more natural to work with the individual chemical 
potentials. Then the final equations retain the natural symmetry between different chemical constituents. 

The ramifications of a lack of partial pressures can be extracted by returning to Eq. (14), and inserting Eq. (39). 
We then get 



_ (40) 

X \ X / 

where we have defined 

7? = /r + nxV,p^ + niViP^ (41) 

and used the fact that 

Vi* = -E("-^*Po+^iViP,^) . (42) 

X 

Letting T^^ i be a sum of a conservative and dissipative piece, i.e. 

T^^i = C-i^ + U^U , (43) 

the consistency relation Eq. (14) becomes 



X \ X / X 



(44) 
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The key point is that we have now isolated terms in such a way that consistency of Eqs. (12) and (13) is guaranteed 
by taking 

fi = Y,fl^ Dh = Y,D^U, C^^i=nip^. (45) 

X X 

Here we recall the discussion of Sec. II D in which it was demonstrated that the total stress tensor T^i is symmetric. 
Since the conservative piece C^i can be shown to be symmetric, by using Eqs. (29) and (36) in Eq. (38), this means 
that we must require D^i to be symmetric. However, as we already stated following Eq. (26), this does not mean 
that we should assume that the individual D^^i are symmetric. 

Finally, the equation of motion for the x*'^-constituent takes the form 

9i< + V, «pf + D^^,) - {n^^^pl + niV^p^) = /f , (46) 

while the total external energy creation rate per unit volume becomes, cf. Eq. (33), 

e-t = ^ (</f + D-UVjvi - [pS + vipf] Fx) . (47) 

X 

One can verify that our equations are consistent with those of earlier formulations, with the exception that our 
formulation has led to a balance of forces, Eq. (45), that docs not include the various dissipation stress-tensors D^^i. 
Thus, to complete the system we can specify each D^^i, and each /f independently, with only the latter having to 
add to the total external force density fi acting on the whole system. We will now show how to limit further the 
possible forms for the /f and D^^i, by employing Onsager's formulation [1] for multi- fluid systems. 



III. THE ONSAGER FORMULATION FOR DISSIPATIVE MULTI-FLUIDS 



Because we have many independent fluids, the number of potential dissipation coefiicients can be quite large. 
However, Onsager [1] demonstrated long ago that microscopic reversibility implies certain equalities among "off- 
diagonal" pieces of the entropy creation rate. Since this argument plays a key role in our analysis it is worth 
outlining the main ideas. A more detailed, pedagogical, description can be found in [44] . 

We begin by noting that the entropy, here represented by the number density rig = s, is maximal for a system 
in equilibrium. This means that any perturbation away from the equilibrium must be represented by quadratic 
deviations. Specifically, in the "thermal frame" associated with the entropy velocity vl, the conservation law Eq. (6) 
implies 



A 



s 



(48) 



At At 

Comparing this to the anticipated expansion, see Jaynes [45] for an elegant exposition of this, near an equilibrium 

s^s^^-^Y.XaL-OXt,, (49) 

a, 6 

we can identify [46] 

1 ^ 

rr, = - - ^ x^L-'Xb = J"^a , (50) 

a,b a=l 

where the N individual Xa are known as "thermodynamic forces" and the J" as "fluxes." The thermodynamic forces 
represent a measure of the departure from global equilibrium in the system, with the fluxes arising in response. The 
Onsager symmetry principle simply states that L°-^ = L^°'. 

Let us now see how we can apply this idea to our formalism. Our entropy creation rate is obtained from Eq. (47), 
by solving for Fg. Noting that the chemical potentials fj,^, which are obtained from 



9nx 



(51) 
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are related to the ji^ via 

Jt = tv'- m^'vl + vipt , (52) 
and that the temperature is given by /x^ = T, we can write 



^ext 



fivl + D\V,vl + rr, + ^ {V^^i- + + D^i.V^wQ . (53) 



XT^S 



We recall that fi and D^i are the total force and dissipation tensors, respectively. We have assumed the locally 
isotropic, manifestly Gahlean-invariant form for £ (i.e. Eq. (36)), and we have defined 

n^ff-\n^T^g,,{vi + vi) . (54) 

We will take the total system to be closed, which means that e'^'^* = and fi = 0. In order to apply (50) we also 
need to work in the entropy frame, which means all velocities will appear as the difl[erence w^s = ^^x ~ ^s- 

The "thermodynamic forces," which drive the system to equilibrium in the sense of Eq. (49), are seen from Eq. (53) 
to be (for x ^ s) 

r {X- = Ai-}, a = l 
Xa=l {X^ = <J, a = 2 , (55) 
{{Xii = V'wl}, a = 3 

implying that the corresponding "fluxes" are (again for x 7^ s) 

{{ Jx = -Fx} , a = 1 

{Jf^-/>}, a = 2 . (56) 

{Jf. = -D\j}, a = 3 

The Onsager formulation for our system is thus 

Jt^ -Y. + ^'?k< + ■ (57) 

The tensorial aspects of the "L" and "L" coefficients are handled by assuming that they can only be constructed 
from combinations of the thermodynamic forces, jj^, w^^^ and V^w^g, and the background geometry terms, i.e. the 
metric g''^ and the volume form Cijk = ^/giijk], where [ijk] is completely antisymmetric with [123] = 1. Moreover, 
in keeping with the spirit of the Onsager expansion and due to the fact that we are supposedly close to equilibrium, 
we consider only those terms that lead to quadratic combinations of the thermodynamic forces in Fg. With these 
restrictions it is clear one would not expect coupling between forces with diflferent tensorial nature. We readily flnd 
that Z/^^ = L^^ = and that we can rule out the use of w^^ and V'to^g in any of the coefficients. The most general 
coefficients within our assumptions that can be written are 

Lxy = 7xy = 7yx , 

^"ijki = C^9ij9ki + rf^ (^ikQii + 9ii9jk - ^5ii5fe;^ + ^(^"^^ <^ijme"' kl 

= C^9ii9ki + rf^"^ {gikQji + 9il93k - '^9ij9ki^ + ^(^^"^ ^ijmi^ ki ■ (58) 

Clearly, each coefficient inherits for its spatial indices the symmetries of either gij or Cijk- There has also been a 
lot of thought put into the coefficient i^J^,;. The key point is that the most general four-index object that can be 
constructed is a linear combination of the three terms gijQku 9ik9jh and gugjk (noting that eijm^'^ki = 9ik9ji - 9u9jk)- 
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Our discussion will be greatly facilitated by introducing the trace-free shear 9*{, vorticity ^^'^ expansion ©xg: 



4s 



We will also need the vorticity vector that is dual to the vorticity, i.e 



1 

2!' 



This allows us to utilize the well-known result that 



1 



(59) 



(60) 



(61) 



Here we should emphasize that the proper definition of vorticity, for which one recovers the Kelvin- Hclmhohz con- 
servation theorem, is in terms of the momentum [2, 20]. In the absence of entrainment the quantity we use and the 
conserved vorticity are identical; but when entrainment is present, the difference is crucial. In terms of the analysis 
presented here, our introduction of the "vorticity" is purely a matter of convenience, as it helps in separating purely 
symmetric (in the spatial indices) from antisymmetric objects. 
The fluxes are now of the form 



Tx = -$](7xyM'' + T'^''eys) , 

/r = -25i,X(^''< + ^'''^ys) ' 
yT^s 

and we see that the entropy creation rate (in the thermal frame) is 



(62) 

(63) 
(64) 



(65) 



It is easy to show, given the explicit transformations in [2] that the chemical potentials are Galilean invariant. Since 
all other terms in the right-hand side of Eq. (65) depend only on velocity differences we see that the entire expression 
is invariant. This is, of course, as it should be. After all, the difference between reversible (TTg = 0) and irreversible 
(TFg > 0) processes cannot depend on the observer. 

There are several different strategies to adopt in demonstrating the positive-definite nature of Fg. Common to all 
is the simple notion of producing terms that look like V'-'" MV, where M is a x symmetric matrix, and y is a 
d-dimensional vector. At that point basic theorems from linear algebra (cf. [47]) can be employed that will yield the 
constraints under which the matrix coefficients lead to positive-definite quantities. In our case, the most natural way 
to gather together terms is to introduce two new vectors (in this linear algebra sense): 



{Oxs} 



{<s}' 

{WL} 



whose dimension is d = 2(A'' — 1). Likewise, the matrices that define the quadratic terms are 



Ml 



{7xy} {r^^} 

|^xy| |^xy| 



(66) 



(67) 



The term containing B^^^By' is obviously already of the appropriate form. Hence, once the necessary constraints 
from linear algebra are obtained on the matrix coefficients, i.e. the dissipation coefficients 7^^^, , etc., Fg > can 
be guaranteed in the entropy frame. 

Finally, we return to the inertial frame, in which the entropy velocity can appear independently of the other 
velocities. If we consider the system to be closed, then the entropy creation rate just picks up the additional term 



12 



proportional to D^i. As D^ij still remains undefined, we can just as well work with £)y . Because Dij is contracted 
with V^v^, and it must be symmetric, the only proposal consistent with our main assumptions is 

Dij = - {C°'Qsgij + 2rf°'gik9jiQf) > (68) 

where Qs and O*^ are defined as in Eq. (59) but replacing w]^^ with v].. Using the previous, linear algebra-based 
line of reasoning, we can determine, in principle, the constraints on the dissipation coefficients and rf°^ that will 
guarantee the positivity of DijV^vi, without spoiling the positive-definite nature of the other terms. 

IV. THE NEUTRON, PROTON, ELECTRON AND ENTROPY SYSTEM 
A. The problem of neutron star dissipation 

At zero temperature, the outer core of a neutron star consists of three interpenetrating fluids: superfluid neutrons, 
superconducting protons, and a highly degenerate gas of normal fluid electrons. Because of electromagnetic coupling 
the charged components lock together on a timescale that is much shorter than, eg., the timescale of stellar rotation 
or oscillation. Hence, the problem usually reduces to a two-fluid system. At a finite temperature, the situation is 
considerably more complex. Most obviously, we need to account also for the entropy. Provided that the system is 
far below the various superfluid transition temperatures, the entropy is associated with the electrons. However, there 
will also be contributions from excitations of the neutron and proton condensates. Accounting for these requires 
more thought. Based on the studies of the analogous problem for superfluid Helium there would seem to be (at least) 
two possible strategics. The difference between them is nicely described by Geurst [9]. 

In the first approach, championed by Carter and his collaborators (see for example [2]), each superfluid component 
is made up of the condensate and a massless gas of excitations. In our formalism this would amount to associating 
an entropy component to both the neutrons and the protons. It is then not difficult to show that the "normal fluid 
density" associated with the the x*'^-constituent is directly proportional to the entrainment coefficient a^'^ (cf. [2]). In 
principle, both the neutrons and the protons will thus contribute to the normal fluid. This effect becomes particularly 
important near the superfluid transition temperatures. Far below the transition temperature, the temperature is 
essentially zero and one would expect both a"** and a^** to be eflPectively zero. Anyway, in order to account for a 
"normal" part of each superfluid we need to begin with a four-fluid system. Otherwise there would be no explicit 
distinction between the entrainment coefficients, eg. between those that couple the fluids to the entropy, and those 
(i.e. q;"p, a"'', and a^'^) that couple the species of particles to each other. As the electrons are not superconducting, 
we can assume that — vl which means that we can set q:°^ = 0. 

The second possibility is perhaps closer to the orthodox approach to superfluids [9, 26], in which one introduces an 
ad hoc separation of the mass density into a "superfluid" density and a "normal" fluid density. In our framework, 
this would involve dividing each density into a piece associated with the condensate and a piece corresponding 
to the quasiparticlc excitations. This philosophy was used in the recent calculation of entrainment parameters by 
Gusakov and Haensel [48]. 

At the end of the day, the separation of a fluid constituent into two pieces is purely formal. One could not 
conceivably separate the superfluid condensate from the excitations in a real physical system. The mathematical 
formalism that we have developed should be flexible enough to allow us to represent both approaches to the finite 
temperature problem. As we will discuss elsewhere, the main issue concerns the interpretation of the various variables 
(especially the entrainment coefficients). 

It should also be noted that, despite neutron stars being self-gravitating bodies, we will not couple gravity to our 
model at this point. In the Newtonian context this is easily done by introducing an external force [2]. We do not 
include this force here as our main focus is on the fluid aspects of the problem. 

Let us proceed by assuming that there are four independent densities, eg. the number densities for the neutrons, 
protons, and electrons, as well as the entropy density. These are denoted (rin, np, no, s). A priori, the system will 
consist of four independent velocities, eg. the neutron, proton, electron and entropy velocities which are denoted 
{vl^,Vp,vl,vl). The entropy fluid is massless. We consider the four-fluid system first, in order to identify all the 
entrainment coefficients, before reducing to the case where all "normal" fluids flow together. We believe that this 
strategy is more natural than one which imposes the existence of only two distinct transport velocities from the 
beginning. 

For the four-fluid system, the "first law of thermodynamics" is given by the expression 

d£ = ^indrin + /UpdUp + /Xedrie + Tds + a^'^dw^p + a^'^dw'^^ 

+ a^Pdw^p + a"Mw;^, + a^^^dwl^ + a'^^dwi , (69) 
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where we are again assuming the Galilean-invariant, locally isotropic form for £. Because of overall charge neutrality 
we have 

rip = rie . (70) 

This, in conjunction with Eqs. (6) and (7), implies 

Tp = Te = -r„ (71) 

and 

=m^ + nf = m , (72) 

where m", wP, and nf are the neutron, proton and electron masses, respectively. The electrons should also closely 
track the protons because of the electromagnetic attraction [17, 27]. Thus, it is a good approximation to simply set 

vi = vi = vi. (73) 

In a later work, we will consider electromagnetic coupling in more detail in order to justify fully this constraint (see 
also [49]). 

In the absence of thermal conduction vl = v].. Because of this the entropy flux is often written as 

svi = svi + ^q' ^ q' = sTwl , (74) 

where 5* is the heat-flux vector. However, by introducing the heat-flux vector we break the natural symmetry 
of our formalism since is put on a special footing. In general is dynamically independent of the other two 
velocities. A standard result is to argue that any deviations of vl from co- motion with will be due to gradients 
of the temperature. Because of the deflnition of temperature, cf. Eq. (69), the entropy flux still represents a 
dynamically independent variable, since T depends on the relative velocities w^yi which include vl- This kind of 
velocity dependence is manifested in the momenta pf, for instance, via the entrainment coefficients a^^ . 
It is convenient to define 

a""" = a"P + a'"' , a"" = a^^ + a"' , (75) 

and 

= + . (76) 
Now the three independent fluid momentum densities can be written 

< = Qij {mn^vi - 2 [(a- + «-) < - «-<] ) , (77) 

< = nppf+n^pt=g^J{mncvi + 2[a'''wi^+a'''wQ) , (78) 
TTf = 2gij (a-«;i, - [a^ + a^] <) . (79) 

Defining 

D'^^i = D^^i + D^^i , (80) 

we find for the three force densities, cf. (46), 



= + V, {vinf + D^h) + ncV, f^c - \mvi\ + tt^V.z;^ , (82) 



(81) 



2 

ft = dtwl + Vj {viwt + D'^h) + sViT + n'^Vivi . (83) 

We can put the finishing touches on our "neutron star model" , by building the thermodynamic "forces" and "fluxes." 
Key to this will be the assertion that the thermodynamic forces are to be considered as linearly independent, i.e. if a 
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linear combination of thermodynamic forces is to be equal to zero, then the coefficients multiplying each force must 
vanish. This assumption does not seem unreasonable. We must first consider the constraints that result from the 
equalities among the particle creation rates Fn, Fp, and Fg. 



By setting Fe 



-Fn, and imposing linear independence, we determine that 

nn — fpn — Jen — /pp — /ep — /ee 



and 



rPP 



np 



If we define 



then the fluid forces are 



where X = 1^ ,A ,Q ,r} ov <j 



fi = - 

and using the fact that we are considering a closed system, we also have 



1 



ft = n = 

Finally, the independent dissipation coefficients become, cf. Eq. (64) 



We also get 



(84) 



(85) 

(86) 
(87) 



(88) 



(89) 
(90) 



(91) 



(92) 

(93) 



(94) 



It is worth noting that this construction ensures that Dij is symmetric, as required for a closed system. 

The final step is to use the positivity of the entropy to set constraints on the dissipation coefficients. It follows 
from Eqs. (65) and (68), and some algebra, that the entropy creation rate is 



TF3 = C°'Ql + '^'n"''9ikgoiQi'Qf 
+2gik9ji [eg es: 



Qkl 
^ns 
Qkl 



+ ens ees] 

+2gij [wg 1^4] 



/nn 


_^nn 


2t""" 






_^nn 




^nc 




e„s 


2t'"" 


^nc 


^cc 







-^nn 


^nc 


A''"" 












A'"= 


A"""" 






A""" 


A""" 


^nn 








A""" 


A"" 











(95) 



This expression shows that the problem is extremely rich. Our system has in general 19 independent dissipation 
coefficients. This can be compared to the general case for superfluid Helium which, according to Putterman [39], 
requires 13 independent coefficients. 
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In principle, it is now a straightforward matter to determine the constraints on the dissipation coefficients, since 
each defining matrix is symmetric and real, and thus linear- algebra theorems can be employed to determine when 
each contribution is positive definite, specifically, if all the eigenvalues of the defining matrix are positive then the 
quadratic form will be positive definite. For example, it is not difficult to show that if both r;"" and 77'^'^ are negative 
then there are no values for the coefficients that will lead to a positive-definite form. On the other hand, if 77"" > 
and 1]'^'' >0, then the quadratic form is positive if and only if 

A similar analysis for the second and third matrices would be tractable analytically, since the characteristic equations 
for the eigenvalues will be third and fourth order, respectively. However, it is not clear that working out the detailed 
conditions would add any further insight at the present time. After all, many of the coefficients in our general 
expression are new, and it may be more productive to begin by seeing if we can interpret their physical meaning. 



B. Final reduction: The two-fluid equations 



In order to connect with previous work on superfluid neutron stars we conclude our discussion by making one 
further reduction. We will neglect heat conduction and let the entropy flow with the charged fluid components. This 
simply means that we take u| = in all equations in the previous section. This implies that we only have two 
independent velocities in the problem, and hence only need two equations of motion. Given this, it is natural to 
combine the momenta of the entropy and the charged components and consider 

TT^i = ttI + < = gij [mricvi + 2 (a""^ + a"') <J = mricgij [vl + ec<c) (97) 

where we have defined the entrainment parameter 

^^^2Kl±i^_ (98) 
mric 

Next, combining Eqs. (82) and (83) we see that 

= -/r-imr„<^-TOrn<-v, (!)=■', + . (99) 

Rewriting this equation in terms of the explicit transport velocities, as in [2, 50], we obtain 

mnc {dt + £vj [v- + ecwf] + n^Vi ^/Xc - ^mv"^^ + sV^T 

= -/f -\rn{l- 2e,) r^M;^ + V,- {D^h - D^) , (100) 

where wc have used the fact that the Lie derivative (which measures the rate at which a quantity changes relative to 
its motion) of a covector follows from 

£yai = v^Vjai + ajViV^ . (101) 

We have also made use of Eq. (6). 

The corresponding equation for the neutrons follows from Eq. (81). We get 



mrin {dt + £vj [vf + SnWr] + ^^nVj ^/in - ^'^^'n^ 



/f + -m (1 - 2£„) rn<" - V,I?"^i , (102) 



where 



16 



The two equations of motion, Eqs. (100) and (102), can be directly compared to Eqs. (176) and (177) in [2]. If we 

introduce 

/r* = t + lrn{l- 2£n) rn<" (104) 

we find that the two results agree perfectly in the zero- temperature (T = 0), non-dissipative {D'^ij = Dij = 0) 



limit. Thus our representation reduces to the familiar result in the relevant limit. In the more general case that we 
have considered, several new aspects appear naturally. Finite temperature effects are explicitly represented via the 
presence of the sVjT term in the equation for the charged fluid/entropy equation. Implicitly, we are also accounting 
for the possibility of a normal part of the neutrons via the entrainment coefficient a"®. The two fluids are coupled 
in a number of ways. In particular, the mutual coupling via /" and -D"y , which accords with Newton's Third Law, 
should be noted. The presence of the "total" dissipation in terms of Dij is also relevant. 

It is instructive to consider the above equations of motion in the context of what we already know about dissipative 
fluids. In particular, we would like to understand which of our many coefflcients can be considered as known, and 
which represent potentially new aspects of the problem. To carry out this exercise, let us begin on familiar territory 
and consider the total dissipation part. From our definitions, we have 

A, = {C°'0s9^J + 2if°'g,k9,iQf) (105) 
which should be compared to the standard right-hand side of the Navier-Stokes equations [26] 

+ V, (CV,i;') = (2r;e,,) + g^,W (CG) . (106) 

Comparing the two expressions we can readily identify the coefficients of the standard shear and bulk viscosities as 

n = r/*°* and C = C*°* • (107) 

Not surprisingly, the interior coupling terms arc not as straightforward to explain. This is quite obvious since we 
have tried to write down a more or less general expression for the permissible terms in a multi-fluid context where 
previous studies have been limited. The problem is exacerbated by the fact that the bulk of the extant work was 
carried out in the "orthodox" framework [26]. To interpret the result, we consider the combination 

+V^' {gij [r"" (Mn - Mc) - C^Onc] + 2r/""e^; + <7"X'} ' (108) 

which follows from the various definitions and some straightforward algebra. Given this expression we can interpret 
the individual terms in the following way: Since it is proportional to the relative velocity, the term involving TZ^^ is 
associated with a force analogous to standard resistivity. Meanwhile, relates to a Magnus-type force which acts 
orthogonally to the relative flow [20]. The two terms involving and Ofj are simply analogues of the standard 
bulk- and shear- viscosities, although here it is the expansion and shear of the relative flow which is important. The 
coefficient t"" accounts for effects due to the fluids being driven out of chemical equilibrium by the flow. In recent 
work Carter and collaborators have refered to this type of interaction as "transfusion". The flnal term in Eq. (108), 
related to cr"'^, appears to be new. 

To conclude our discussion, let us consider the particular form of the mutual friction force which arises because of 
a balance between i) the Magnus force acting on quantised superfluid neutron vortices and ii) scattering of electrons 
off of the vortices. As we have shown elsewhere [20] , this results in a force 

= B'p.Mjkn^w':^ + Bpnn,e,,ke'""'k^Kiw^: , (109) 

acting on the neutrons. Here riy is the vortex surface density and is a vector which is aligned with the vortices 
and has magnitude h/2mn (a hat denotes a unit vector); see [20] for further details. 

Comparing this expression to Eq. (108) we see that only the first two terms in the latter are relevant. In order to 
incorporate Eq. (109) in our framework, we need to be able to identify the coefficients 7?."" and Vi^"". We easily 
see that we must have 

Vi^"" = -B'p^n^Ki . (110) 
The corresponding expression for TZ™ follows from 

fip„n„e,,fee'=''"«^K;<" = Bp^n^n {k,kj - gij) w^^ = 27e'^"<" . (Ill) 
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Multiplying each side by wl^ we find that 




(112) 



The above expressions can now be used, given the explicit estimates for the two coefficients B and B' obtained in 
[20] , to incorporate mutual friction in our dissipative neutron star model. This will allow us to extend our previous 
studies of superfluid neutron star oscillations and instabilities in a number of exciting directions. 



We have described a general formaUsm which can be used to model dissipative multi-fluid systems. Our discussion 

was based on the key conservation laws for mass, energy and Hncar and angular momentum, together with recent 
progress in deriving the equations of motion for the non-dissipative problem from a variational principle [2-5], and 
the classic Onsager symmetry principle which provides a strategy for including dissipative terms in such a way that 
the second law of thermodynamics is satisfied. Following this approach, we derived a set of multi-fluid equations 
which has a number of advantages over the "standard" formalism used to model, in particular, superfluid He^. The 
most important difference is that we have made due distinction between velocities and momenta, which means that 
the entrainment effect is accounted for in a natural way. 

As an application of the general formalism, we considered two models for a neutron star core. The first model 
accounts for three distinct fluids, the superfluid neutrons, a conglomerate of charged particles (protons and electrons) 
and the entropy. We have shown that this model, in principle, requires 19 distinct dissipation coefficients to be 
determined. Detailed microscopic analysis is needed to understand the relevance of the majority of these coefficients. 
We also discussed how the normal fluid fraction of neutrons that should be present at finite temperatures can be 
associated with the entrainment parameter a^^ which encodes how the equation of state depends on the relative 
motion between the neutrons and the entropy. Our second model is obtained by neglecting heat conductivity. This 
couples the entropy to the charged fluid and reduces the number of independent degrees of freedom to two. We 
discussed the form of the various dissipative terms in this two-fluid model, and interpreted their physical meaning. 
We also showed how the mutual friction, which is mediated by the superfluid vortices, can be accounted for within 
our model. 

One can envisage a number of interesting applications of the framework that we have developed in this paper. In 
the case of neutron stars, the possibilities range from the viscous damping of neutron star oscillations due to electron- 
electron scattering [19] to intricate issues concerning the effect of entrainment on the mutual friction force [20]. Since 
our model allows for much more complicated multi-fluid dynamics than has so far been considered in the literature, one 
woiild have to be somewhat suspicious about any statements about the damping of superfluid neutron star oscillations 
which were made without consideration of the various degrees of freedom. It would certainly be interesting to revisit 
the problem of oscillation modes driven unstable by gravitational radiation [11] and investigate the role of the true 
multi-fluid degrees of freedom. Our hope is this work will inspire a significant improvement of our understanding 
of superfluid neutron star dynamics, and perhaps generate input from other communities (eg. chemistry) that have 
much experience with multi-fluid systems. 

Finally, it is worth emphasising that our formalism is quite general, and one would expect it to be useful also 
in other problem areas. Applying it to other multi- phase flow problems, as discussed in for example [6-8], seems 
relatively straightforward. There are also more exotic possibilities. In particular, we believe that our formalism could 
be quite relevant for problems involving fluid analogues of various curved spacetimes, such as black holes and big-bang 
cosmologies, that can be potentially observable in laboratory systems (see [51] for a discussion of "sonic" horizons in 
superfluids, and [52] for a recent review). For this kind of application, it would be particiflarly important to establish 
whether the "spacetime" effects are likely to be rendered experimentally undetectable by various forms of dissipation. 



V. 



CONCLUDING REMARKS 
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Appendix A 



Here we will motivate why the conservative part of the total stress tensor C^i can be obtained from a variation of 
the action with respect to the metric. The analysis is by no means rigorous, but it docs lead to the result that has 
been obtained rigorously in Ref. [2] . Our aim is to give the reader not familiar with all aspects of action-based field 
theory some basic physical intuition, using simple energy conservation arguments. 

Let us consider a displacement in a small amount of time 5t = t2 — ti of a closed, isolated system such that the 
individual constituent particle numbers as well as their corresponding velocities are kept fixed. We will also suppose 
that each constituent undergoes the same displacement so that their respective velocities can be written 

St 

Because the system is closed and isolated, the total energy will be conserved. Prom Eqs. (15) and (28) the change in 
total energy can be written 

5U = U{t2)-U{h) 




^ -S [ CdV, (114) 

Jv 



dt 



where the last term represents the difference in the integral before (at time ti) and after (at time t2) the displacement 
is put in place. If we now impose that the total energy does not change 6U = 0, and that the velocities remain fixed, 
i.e. dvl^/dt = 0, then upon insertion of Eq. (13) we find 

5 1 £dv = J2 r [ <(7?-v,c"^^0 dV 

Jv ^ Jti Jv 

{h - ^jCh) dV , (115) 



where the second equality follows because the displacement is the same for each constituent, we have used Eq. (14), 
and we have kept only terms linear in (5^*. Since the system is isolated, the net force = 0. Finally, after integrating 
by parts, we can write 

5 ( CdV^l [ C'^ {WiSij + Vj<56) dy - / ChrijSC dA . (116) 

Jv ^ Jv JdV 

The last term will vanish if we impose suitable boundary conditions on the control volume (or assume that the 
boundary is well outside the fluids). 

We now need to establish that the deformation of the system set up by the displacement (5^* induces a change in 
the metric, which in fact will define the 6 variation in Eri. (116). Consider an active coordinate transformation where 
the metric is pushed (at time to the new points = x* + J^*. Denoting the transformed metric as g^j, then 
ordinary tensor analysis implies 



_ _ dx'^ dfx'' 

Next we define 



^9^1=9^Ji^'')-9^Ax'"), (US) 

and, keeping terms up to linear order in the displacement, we find 

Sgij=S/id^j + Wj5^i . (119) 
Thus in Eq. (116) we can rewrite the displacement on the right-hand side in terms of dgij. 
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which implies 



Next we must determine how to vary the fundamental variables and under the displacement. Recall that we 
are assuming that the number of x"^-constituent particles A'x remains fixed. This implies 

6N^ = 6 [ rix dy = Sn^ = -n^^ , (120) 
Jv Vd 

where matrix theory (cf. [47]) can be used to show 

SV9=lV99''Sgij . (121) 

Thus, using also Eq. (119), we know how the displacement changes rix. Furthermore, we also know how the particle 
number density currents change since = n^vl., but the velocities remain fixed, and so ^n^ = i;*5nx. Therefore, 

we have 

Sn^ = —n^g'^Sgij , (5< = ~Kg^^5gjk ■ (122) 

The small variation indicated on the left-hand side of Eq. (116) is the difference between the integral before the 
displacement and the integral after the displacement. But as the integral is a function only of rix, n]^, and gij, and 
the variations Jrix and dnt^ induced by the displacement are equivalent to a variation of the metric, we see that the 
5 on the left-hand side of Eq. (116) can be equally understood as a variation of the metric. Thus, Eq. (116) can be 
reduced to 

'c--A^),,,dy = 0, (123) 
V \ Va 9gij J 

= (124) 

This result should be compared to the standard formula used in general relativity (see, for example, [42]). 
Finally, using the less-general internal energy given in Eq. (36) we find 

(7^^ = vl/^^^+^nip^ (125) 

X 

which is Eq. (38) in the main body of the paper. 

Appendix B 

Consider the simplest application, where each particle number is conserved independently and there are no forces 
acting on a constituent, i.e. /f = 0. Then 

dtn^ + Vi< = (126) 

and 

I M) + (p^ni) = nxViP^ + niViP^ . (127) 

Recall that 

Vi<i> = -J2i^-'^iPo+KVip]) . (128) 

X 

We see that the term appearing in the summation is precisely the term on the right-hand-side of the above "Euler" 
equation. Hence, unless ^ takes the particular form of 

JV 

* = ^\plpl)+i'\plpl) + ... + V^(Po^,pf ) = EV'"(fo,K) , (129) 

CT=1 
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which results in the exact differential 

= - {Ua^iPo + K\/ip'^) (130) 

for each constituent, such that 

dip" i dip" 



(131) 



the right-hand-side of the "Euler" equation will not be an exact differential, and so cannot be written as a total 
divergence. 
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